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GRID GENERATION FOR TURBOMACHINERY P ROBLEMS 


1. Introduction 

In this report we outline the development of a computer code to generate numerical 
grids for complex internal flow systems such as the fluid flow inside the space shuttle 
main engine. For computing the viscous compressible flow inside the main engine where 
the fluid undergoes complex turns through various ducts, it is necessary to discretize the 
physical space while fulfilling several competing requirements. The topology of the grid 
structure should depend on the flow solver algorithm being used. Finite difference codes 
usually require constraints on the grid structure such as separability of the indices for 
efficient computational procedures, while finite element codes can be implemented with 
less stringent requirements on the index structure. The grid should provide reasonable 
resolution of the flow field relative to the total number of nodes being used. This requires 
providing more grid points in regions of large gradients of flow qualities, such as the 
viscous zones near solid boundaries. The numerical grid should meet certain smoothness 
requirements so that the metrics of the curvilinear grid can be computed numerically and 
these computed metrics are continuous. Unreasonably skew’ed grid cells and singular points 
in the grid where the local transformation of the physical space to computational space 
has very small or very large Jacobians should be avoided if at all possible. Otherwise, such 
grids will require special handling by the flow solver algorithm and also may give rise to 
numerical inaccuracies and instabilities. 

There are several grid generation techniques and special purpose codes which can gen- 
erate reasonable grids for simple two dimensional and three dimensional geometries, for 
both internal and external flows. These techniques fall under two classes: algebraic gener- 
ators and elliptic generators. .Algebraic generators use various interpolation and stretching 
functions while elliptic generators solve a set of elliptic partial differential equations. While 
both techniques are effective for simple geometric regions, it is usually difficult to use them 
to develop a composite grid over a complex internal flow domain. It is particularly dif- 
ficult to enforce a smooth transition from one zone to another if the physical domain is 
composed of separate pieces with different geometries. In this report we describe a new 
blending method to generate smooth grids over complex domains by systematically blend- 
ing several smooth grids, each of which is developed to conform to different pieces of the 
domain. In particular, the development of a grid for a 3-D rotor-stator combination in 
turbomachinery will be outlined. 

The first step in developing a grid is to decide upon a grid index structure, which 
amounts to choosing a global strategy for mapping the physical domain boundaries to the 
boundaries of the computational domain. The total domain is broken into several simpler 
overlapping domains. Each subdomain should have a part of its boundary in common with 
a part of the physical boundary of the total domain. The rest of the boundaries of the 
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subdomain are not restricted and are chosen according to convenience. A suitable algebraic 
grid is generated on a subdomain which conforms to the grid spacing and smoothness 
requirements of the local physical domain boundary. When all the subgrids are generated 
by simple methods, they will overlap each other and conform to different parts of the 
physical boundary of the total domain. Suitable weighting functions are then developed 
to blend all the grids together to generate a smooth grid over the total domain, which 
conforms to the various boundary requirements. This technique has been developed and 
tested by Steinhoff for complex two dimensional^^^ and three dimensional^^^ geometries. 

In this report we outline the blending technique for generating a grid for stator- 
rotor combination at a particular radial section. The computer programs which generate 
these grids are listed in the Appendices. These codes are capable of generating grids at 
different cross sections and thus providing three dimensional stator-rotor grids for the 
turbomachinery of the space shuttle main engine. 

2. Blending Methods for Grid Generation 

In many cases where a smooth computational grid is required, the boundary of the 
computational domain can be decomposed into a number of pieces, each of which is fairly 
simple. We suppose that an adequate grid can be easily generated for each of these pieces, 
if considered by itself, and describe a method for blending these "elementary'’ grids into 
one smooth composite grid which has all of the pieces as its boundary. Examples where this 
technique can be used include external flow over a entire aircraft, where simple methods 
exist for generating grids individually over each of the lifting surfaces and the pieces of 
the body. Other examples include internal flows in turbomachinery where methods exist 
for generating grids for each element such eis a rotor or a stator taken separately. ,\n 
important feature of the concept is that it can be used recursively. Composite subgrids 
can first be formed from elementary grids, using the method, then, the same method can be 
used to form larger composite grids out of these individual subgrids. If algebraic methods 
are used to form each elementary grid, which can often be done since each piece is simple, 
then the entire grid generation procedure is algebraic, since the blending is non-iterative 
and involves no partial differential equation solutions, .\ccordingly. where applicable, it 
is a fast method suitable for interactive use. Also, if a partial differential equation is to 
be solved for some physical quantity and an iterative method is used to solve a set of 
discrete equations on the grid, which is usually the case, then at each iteration the grid 
can be quickly re-generated and there is no need to store the entire grid system. This 
feature can be especially important for large three dimensional problems. This method is 
very different from other algebraic methods, such as those of Eisemen^^l. Each elementary 
grid is taken to be previously determined, either by algebraic methods, partial differential 
equation solution or any other means. These grids can be defined over the entire space, 
rather than just on surfaces as in ■‘transfinite interpolation” schemes. 
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An important feature of the method is that it allows the grid designer to use software 
packages and methods already developed or being developed by others (which can be quite 
sophisticated and complex) for the elementary grids about each piece of the problem. 
These can be used as “black boxes,” and after each elementary grid is generated the grid 
designer can blend them together. Also, after a composite, complex grid is generated, if 
one of the pieces is later modified, only the single new elementary grid need be recomputed 
and blending into the composite grid. 

In this report, grid generation for a stator-rotor combination that exists in the space 
shuttle main engine is described. Elementary //-grids are formed for the basic profile 
shapes of the stator and rotor by algebraic methods, they are each blended with proper 
outer boundary grids, each is given its own camber and rotation and then the composite 
grid is obtained. Figures of the grids for one radial cross section are shown, but the final 
computer program will generate three dimensional stator-rotor grids with as many cross 
sections as necessary. 

2.1. The Basic Method 


Consider a set of N grids, each spanning the same computational space and approxi- 
mately the same physical space. For simplicity, we define the computational coordinates to 
be just the (integer) indices of the grids. Thus, in n dimensions we have an n component 
vector, Fm{i){= (jJm(0» l/m(^)i ~m(^)) for n = 3) defined on each grid (labeled m) as a 
function of the indices £(= {i,j,k) for n = 3). It is important to think of the n compo- 
nents of fm as ordinary smooth functions defined in the computational {(.) space. Defining 
non-negative weighting functions P"*(£), the physical coordinates of the composite grid 
are then simple weighted sums of those of the elementary grids: 



Yi P"'WF„{e) 


/ 




The weighting functions are, in general, functions of all of the indices £, and are a func- 
tion of how close the point i is to the elementary surface segments. When t approaches 
some surface segment, say mi, then P'^‘(£) must approach 1 and all the others P’s must 
approach 0 since there we must have 


r^(^) — * 

Some of the “art” of using the method resides in the determination of the functions 
P”*(£). Since values of r,„(£) which define smooth grids are determined separately about 
each elementary surface, the P'”(£) do not have to do as much work as in an interpolation 
method where they typically completely determine one of the coordinates. In the next 
section, it will be seen that very simple functions are sufficient. The main problems arise 
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when grids must be blended with very different values of r in certain regions of £ near 
an elementary surface. Then, care must be taken that a number of derivatives of P'^[t) 
are 0 as £ approaches the elementary surface (mi), in addition to the value of (£) 
approaching 1. As more derivatives are made to go to 0, the region in (£) space where 
fc{t) approaches rmi(^) becomes larger. 

2,2. ff-grid For a Stator- Rotor Combination 

To develop /f-grids for stator-rotor combinations, a sequence of elementary grids are 
generated and blended together sequentially, /f-grids are developed for a stator and a 
rotor separately by the same method. For example, the coordinates of the profile shape 
of a stator cross section are input into a subroutine to fit a cubic function for the mean 
camber line of the profile. The cubic is used to define a vertical shearing to approximately 
straighten the airfoil. After the grid is generated this shearing will be applied in reverse to 
all the grid points so that the initial airfoil is recovered. The shearing function (of j) is a 
straight line in front of the leading edge and behind the trailing edge, matching the slope 
and position of the mean camber line there, and is an interpolating cubic function of i in 
between. This function is simply subtracted from the initial airfoil coordinates and, after 
the mappings are complete, added back to each of the grid points to generate the final grid. 
A ba^ic /f-grid is generated about the profile shape with camber removed by an algebraic 
method developed by Pelz and Steinhoff^°^ A detailed study of this mapping was presented 
in Ref. (5) for a single airfoil, where it was shown that a particular transformation can 
be used to eliminate the singularity which normally arises at the leading edge in this case. 
A compressible flow problem was solved on this grid and the solution was showm to be 
accurate once this singularity was removed. Program I, listed in .Appendix I generates 
basic H-grids by this method. A coarse ff-grid generated by this code is shown in Fig. 
1 for a stator profile with camber removed. This grid is blended with a 2nd grid which 
is a rectangular Cartesian grid with a slit in the middle and has the same topology and 
number of cells as the basic //-grid. The objective is to blend grids 1 and 2 to generate a 
grid that will approach the H-gnd near the airfoil and the rectangular grid near the outer 
boundary. Let (fi,/ 2 ) and ( 71 ,^ 2 ) be the interval of i and j indices representing the airfoil 
surface in the iJ-grid. Let (foifs) and {jo,h) be the intervals of i and j indices for the 
outer boundary. The blending is done with a single weight function. p(£) : 

= P(0n(£) - (I - P(^))r'2(f~) 

where 

p(£) = ~T - cos (jra)! [l - cos (^rd)} 

4 ‘ 

a(i) = (i - to)/(»j - I’o) , I’o < i < I'l 

=1 , *1 < * < *2 

= («3 - ‘)/(»3 - *2) , »2 < ‘ < »3 
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j3 function is defined similarly. The resulting blended grid is shown in Fig. 2. The outer 
boundary of the blended grid has the appropriate spacing such that it can obey periodic 
conditions on the top and bottom and it matches with rotor or stator grids on the sides. 
The inner ^f-grid near the airfoil is good for some flow solvers. However, in order to 
provide better resolution near the leading edge and to simulate the blunt trailing edge, 
the mesh lines emanating from the leading and trailing edges are split into two lines, thus 
adding another row of points in the grid. Now the camber is added to all the points and 
the resulting grid is smoothed with a simple Laplcian type of filter. The computed grid 
(40 X 18) for a stator is shown in Fig. 3. Fig. 4 shows a fine mesh (40 x 34) for the stator. 
For a rotor grid, the coordinates are rotated about the leading edge by an angle (32.5° in 
this case). Figs. 5 and 6 show coarse (40 x 18) and fine (40 x 34) grids for a rotor. 

Finally the origin of the rotor coordinates are shifted and the stator-rotor grids are 
matched to derive a composite grid. Fig. 7 shows the composite stator-rotor grid with 
coarse spacing (80 x 18), repeated for several blades. The spacing of the exit boundary of 
the stator grid matches with the spacing of the inflow boundary of the rotor. Thus the two 
grids always match if they slip past each other in whole increments of the mesh spacing. 
This feature is necessary for some flow solver algorithms. 

Appendix II contains the listing of Program II which takes the output of Program I, 
namely the basic ff-grid for a stator or a rotor profile with camber removed and produces 
the final stator-rotor grids described above. 
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Figi, 1 Basic H - Grid for a Stator 
Profile with Camber Removed 
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H - Grid blended with an outer 
rectangular grid 
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Rotor Grid - Fine Mesh (40 
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Fig. 7 Composite Stator - Rotor Grid (80 





Appendix I 


CCC 


C 


C 


PARAMETER (kx=189,LY=34) 

IMPLICIT REALJkB (A-H,0-2) 

COMMON XU(KX,LY),XL(KX,LY),YU(KX,LY),YL(KX,LY), 

1 SCAL,NX,NY,ILE,ITE 

COMMON/MSH/ AO(KX) ,B0(&7) ,XSING,YS ING, XL IM, BOUND, HWALL 
COMHON/GEO/ XP(KX) , YP(KX) , TRAIL, SLOP!, NP 
COMMQN/FLO/ FHACH, ALPHA, CA,SA,C IRC, KSYM 
COMMQN/RLX/ ERES , ARES , DG , AG , PI ,P2 ,P3 ,P4 , P5,U IS,QC, 

1 IRES,JRES, IQ, JG,NSUP, KOBE, MODE 

COMMON/OUT/ SVU(&5) ,SVL(65) ,SMU(6S) ,SML(65) ,CPU(65> ,CPL(65) 
DIMENSION T0T0(6),C0V0(6),P10(6),P20(6) ,P30(6),P40(6>,P50(6) , 

2 GMESH(6),TITLE(20) ,RES(502) , COUNT (502) 

IREAD = 5 

lURII = 6 

RAD = 57.2S5779513082 

READ (IREAD, 530) (TITLE( I) , 1=1,20) 

WRITE (IURIT,630) (TITLE( I) , 1=1,20) 

READS PROFILE DATA FOR CAMBER REMOVED CASCADE 
READ (IREAD, 500) 

READ (IREAD, 510) ESYM,FNU,FNL,FNX,FNY,FMESH,VIS,QC 

ISYM = FSYM 

NU = ENU 

NL = ENL 

IF (NU.LT.l) GO TO 301 

NX = ENX 

NY = ENY 

MMESH = EMESH 

IE (QC.LE.O.) QC = .9 

READ (IREAD, 500) 

DO 4 M=l, MMESH 

READ (IREAD, 510) TOTO(M) ,C0V0(M) ,P10(M) ,P20(M) ,P30(M> , 

1 P40(M),P50(M),GMESH(M) 


4 CONTINUE 

READ (IREAD, 500) 

READ (IREAD, 510) FMIT1,EMIT2,H«ALL 
MITl = FMIIl 

MIT2 = FMIT2 

READ (IREAD, 500) 

READ (IREAD, 510) ALl ,AL2,SIEP,FM1 ,FM2,DMACH 
IF (STEP.LE.O.) AL2 = ALl -STEP -1. 

IF (DMACH.LE.O.) EM2 = EMI -DMACH -1. 

AL = ALl 

EMACH = EMI 

READ (IREAD, 500) 

READ (IREAD, 510) TRAIL, SLOPT,XSING,YSING 
NP = NL +NU -1 

READ (IREAD, 500) 

DO 12 I=NL,NP 

READ (IREAD, 510) XP(I),YP(I) 

XP(I) = XP(I)A2.0 
12 CONTINUE 

L = NL +1 

IE (ISYM.GT.O) GO TO 15 
READ (IREAD, 500) 
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oooooorjoo 



DO 14 

1=1, NL 


READ 

(IREAD,510) 


J 

= L -I 


XP(J) 

= VAL 

14 

YP(J) 

= DUM 


GO TO 

21 

15 

J 

= L 


DO 16 

I=NL,NP 


J 

= J -1 


XP(J) 

= XP(I) 

16 

YP(J) 

= -YP(I) 


21 CHORD = XP(1) -XP(NL) 
ALPHA1=0./RAD 
L1=NL+1 
J=L1 


ISYM=0 

XM = XP(NL) +.25ACH0RD 

URITE (IURIT,32) 

32 FORMATdSHOINPUT COORDINATES/ 

1 15H0 X ,15H Y 

DO 34 1=1, NP 

34 URITE (IURIT,610) XP(I),YP(I> 

URITE (IURIT,3&) 

36 F0RMAT(15H0 XSINQ ,15H YSING 

1 15H TE ANGLE ) 

URITE (IURIT,610) XSING, YSING, SLOPT, TRAIL 
TRAIL = TRAIL/RAD 

IF (MHESH.EQ.l) GO TO 51 
DO 42 M=2,MMESH 
NX = NX +NX 

42 NY = NY +NY 

51 CALL COORD 
CALL GEON 


1220 


c 

c 


64 


c 

c 


CALL HESH 

URITE(6,1220) SCAL,ITE,ILE 
URIIE(6,A) CHORD 
F0RMAT(E15.8,2I3) 

URITE (IURIT,600) 

URITE (IURIT,62) 

DO 64 1=1, MX 
XA = SCALAXUd,!) 

YA = SCAL*YU(I,1) 

XB = SCAL*XU(I,JU) 

YB = SCALAYU(I,JU) 

URITE (IURIT,680) I, AO(I) ,XA,YA,XB,YB 
URITE (IURIT,600) 

URITE (IURIT,66) 


66 F0RMAT(5H 
1 
2 
3 


500 


KY 

STOP 

FORMAT(IX) 


J,15H 

15H 

15H 

15H 

NY +2 


BO , 
X<1,J), 
Y(1,J),15H 
Y(MX,J)) 


) 


,15H TE SLOPE , 


X(MX,J), 


15 


510 

F0RMAT<8F10.6) 

530 

F0RMAT(20A4) 

600 

FORMAT(lHl) 

610 

F0RMAT(F12 

.4,8F15.4) 

630 

FORHATdHO 

,20A4) 

640 

FORMATdS, 

7115) 

650 

FORHAK 17, 

E12.4,2I5,2E12.4,2I5,E12.4,2F10.4,I7) 

660 

F0RMAT(E12 

.4,E15.4,3E15.4) 

670 

F0RMAT(F10 

.6,816) 

680 

FORMAT (15 

,7F15.4) 

700 

FORMATdSHOCOMPUTING TIME,F10.3,10H SECONDS) 


END 



SUBROUTINE COORD 


PARAMETER 

(kx=189,LY=34) 

CCC 

IMPLICIT 

REALA8 (A-H,0-Z) 


COMMON 

XU(KX,LY),XL(KX,LY),YU(KX,LY),YL(KX, 


1 

SCAL,NX,NY,ILE,ITE 


COMMON/MSH/ AO(KX) ,B0(67) ,XSING,YSING,XLIM, BOUND 

C 

SET UP THE GRID BOUNDARIES AND DETERMINE NUMBER 

C 

UPSTREAM 

ON THE BODY AND DOWNSTREAM 


AY 

= 1. 


AX 

= 1. 


XLIM 

= 0.5 


BOUND 

= 0.95 


XMAX 

= XLIMABOUND 


SY 

= .5 


YMAX 

= HUALLABOUND 


S 

= 1. 

CCC 

S 

= 0.5 


LX 

= NX/2 +1 


MX 

= NX +1 


DX 

= 2. ABOUND/NX 


L 

= I.OOOOOIAXMAX/DX 


ILE 

= LX -L 


ITE 

= LX +L 


DO 12 1=1 

,MX 


D 

= (I -LX)ADX 

C 

FOR NO STRETCHING COMMENT THE NEXT SIX LINES 

CC 

IF (ABS(D).LE.XMAX) GO TO 12 

CCC 

B 

= 1. 

CCC 

IF (D.LT. 

0.) B = -1. 

CCC 

A 

= 1. -((D -BAXMAX)/(1. -XMAX))AA2 

CCC 

C 

= AAAAX 

CCC 

D 

= BAXMAX +(D -BAXMAX)/C 

12 

A0(I) 

» D 


JU 

= NY/2 +1 


MY 

= NY +2 


DY 

= 2. /NY 


DO 22 J=1 

.,JU 


D 

= (JU -J)ADY -0.5 


SGN 

= SIGN d.,D) 


OF POINTS 


C FOR NO STRETCHING S=1.0 (J-DIRECTION) 

C 22 BO(J) = YMAX*.5*(SGNA(2.ASGN*D)*AS +1.) 
22 BO(J) = YMAXA.5*(SGN*<2.ASGNAD)AAS +1.) 


16 


DO 24 J=1,JU 


J1 = JU 

+J 

J2 = JU 

+ 1 -J 

24 BO(Jl) = -B0(J2) 


RETURN 


END 



c 

c 


SUBROUTINE GEOH 

PARAMETER 

(kx=189,LY=34> 

IMPLICIT : 

REALA8 (A-H,0-Z) 

COMMON 

XU(KX,LY),XL(KX,LY),YU(KX,LY) ,YL(KX,LY), 

1 

SCAL,NX,NY,ILE,ITE 

COMMON/MSH/ AO(KX) ,B0(67) ,XSING,YSING, XLIM, BOUND, HUALL 

COMMON/GEO/ XP(KX) ,YP(KX) , TRAIL, SLOPI ,NP 

COMMON/MPl/ XS(KX) , YS(KX) ,D1 (KX) ,D2(KX) ,D3(KX) 

PI 

= 3.1415926535898 

MX 

= NX +1 

SCAL 

= 2.000002AXLIMAB0UND/(XP(NP) -XSING) 

ANGL 

= ATAN<SLOPT) 

ANGLl 

= ATAN2((YP(1) -YSING),(XP(1) -XSING)) 

ANGL2 

= ATAN2((YP(NP) -YSIHG) , (XP(NP) -XSING)) 

ANGLl 

= ANGL -.5A(ANGL1 -TRAIL) 

ANGL2 

= ANGL -.5A(ANGL2 +TRAIL) 

T1 

= lAN(ANGLl) 

T2 

= TAN(ANGL2) 

ANGL 

= PI +PI 

U 

= 1. 

V 

= 0. 

DO 12 1=1 

,NP 

XA 

= XP(I) -XSING 

YA 

= YP(I) -YSING 

ANGL 

= ANGL +ATAN2((UAYA -VAXA),(UAXA +VAYA)) 

R 

= SCALASQRT(XAAA2 +YAAA2) 

U 

= XA 

V 

= YA 

R 

= SQRT(R) 

XS(I) 

= RAC0S(.5AANGL) 

12 YS(I) 

= RASIN(.5AANGL) 

SCAL 

= l./SCAL 

CALL SPLIE (1,NP,XS,YS,D1,D2,D3,1,T1,1,T2,0,0.,IND) 


RETURN 

END 

C 

C 

SUBROUTINE HESH 
PARAMETER (kx=189,LY=34) 

CCC IMPLICIT REALA8 (A-H,0-Z) 

C GENERATES HESH 

COMMON XU<KX,LY),XL(KX,LY),yU(KX,LY),YL(KX,LY), 

1 SCAL,NX,NY,ILE,ITE 

COMMON/MSH/ AO(KX) ,B0(67) ,XSING, YSING,XLIM, BOUND, HMALL 
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COHHON/GEO/ XP(KX) ,YP(KX) , TRAIL, SLOPT,NP 
COMMON/MPl/ XS(KX) ,YS(KX),D1(KX),D2(KX),D3(KX) 

DIMENSION P(KX) ,Q(KX) ,S(KX) 

PI = 3.1415926535898 

MX = NX +1 

MY = NY +2 

JU = NY/2 +1 

XO = XSING/SCAL 

YO = YSING/SCAL 

1221 F0RMAT(4F10.6) 

XIE = SQRT(XLIM*BOUND +XL1MABOUND) 

YTE = YS(NP) 

XFF = SQRT(AO(MX) +XLIMABOUND) 

A1 = .1 

WRITE(20,1012) MX, MY 

1012 F0RMAT(2I8) 

1013 F0RMAI(2F12.6) 

DO 16 J=1,MY 
II =0 

DO 10 1=1, MX 

A = AO(I) +XLIMABOUND 

R = SQRKAAA +B0( J) ABO( J) ) 

C TYPEA,J,I,R,XIE,A1 

C COMPRESS THE GRID NEAR L.E 

IF(ABS(R).LE. 0.00001) GO TO 1014 
R = ((XTE +A1)/XTE)AR/(SQRT(R) +A1) 

GO TO 1015 

1014 R=0.0 

1015 CONTINUE 

ANGL = 0. 

IF (R.GI.O.) ANGL = .5AATAN2(B0( J) , A) 

IF (J.LE.JU. AND. ANGL. LT.-.25API) ANGL = ANGL +PI 
IF (J. GT.JU. AND. ANGL. LT..25API) ANGL = ANGL +PI 
P(I) = RAC0S<AN6L) 

Q(I) = RASIN(ANGL) 

IF (ABS(P(D).LE.XTE.AND.Il.EQ.O) II = I 
IF (ABS(P(D).LE.XTE) 12 = I 
10 S(I) = 0. 

IF (Il.NE.O) CALL INTPL ( II, I2,P,S, 1 ,NP,XS,YS,D1 ,D2,D3,0) 

DO 16 1=1, MX 

Q(I) = Q(I) +S(I) 

IE (J.GE.JU+1) GO TO 14 

XU<I,J) = XO +P(I)AP(I) -Q(I)AQ(I) 

YU(I,J) = YO +2.AP(I)AQ<I) 

GO TO 16 

14 32 = MY +1 -J 

XL<I,JJ) = XO +P(I)AP(I) -Q(I)AQ<I) 

YL(I,JJ) = YO +2.AP(I)AQ(I) 

16 CONTINUE 
MY2=MY/2 

C OUTPUT GRID COORDINATES FOR USE IN QRID2 

DO 551 3 = 1, MY2 
DO 551 I=1,MX-18 
WRITE(20,1013) XU<I,3),YU(I,J) 
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non 


551 


CONTINUE 
DO 550 0=1, MY2 
DO 550 I=1,MX-1B 
WRITE(20,1013) XL(I,J),YL(I,J) 

550 CONTINUE 
RETURN 
END 

SUBROUTINE SPLIE(M,N,S,F,EP,FPP,FPPP,KM,VM,KN,VN,MODE,FQM, IND) 
CCC IMPLICIT REALA8 (A-H,0-Z) 

SPLINE FIT 

INTEGRAL PLACED IN FPPP IF MODE GREATER THAN 0 
IND SET TO ZERO IE DATA ILLEGAL 
DIMENSION S(l),F(l),EP(l),FPPa),FPPP(l) 

IND = 0 

K = IABS(N -M) 

IF (K -1) 81,81,1 
IK = (N -M)/K 

I = M 

0 = M +K 

DS = S(J> -S(I) 

D = DS 

IF (DS) 11,81,11 
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DF 



= (E(J) -F(I))/DS 


IF 

(KM 

i - 

2) 12,13,14 

12 

U 



= .5 


V 



= 3.A(DF -VM)/DS 


GO 

TO 

25 


13 

U 



= 0. 


V 



= VM 


GO 

TO 

25 


14 

U 



= -1. 


V 



= -DSAVM 


GO 

TO 

25 


21 

I 



= 3 


0 



= 0 +K 


DS 



= S(J) -S(I) 


IF 

(DADS) 

81,81,23 


23 

DF 

= (F(0) -F(I))/DS 


B 

= l./(DS +DS +U) 


U 

= BADS 


V 

= BA(6.ADF -V) 

25 

FP(I) 

= U 


FPP(I) 

= V 


U 

= (2. -U)ADS 


V 

= 6.ADF +DSAV 


IF (0 

-N) 21,31,21 

31 

IF (KN 

-2) 32,33,34 

32 

V 

= (6.AVN -V)/U 


GO TO 

35 

33 

V 

= VN 


GO TO 

35 

34 

V 

= (DSAVN +FPP(I))/(1. 

35 

B 

= V 


= DS 
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noon non 


41 


51 


G1 

71 


DS = S(J) -S(I) 

U = EPP(I) -FP(I)AV 

FPPP(I) = (V -U)/DS 

FPP(I) = U 

FP(I) = (F(J) -F(I))/DS -DS*(V +U 

V = U 

J = I 

I = I -K 

IF (J -H) 41,51,41 

I = N -K 

FPPP<N) = FPPP(I) 

FPP(N) = B 

FP(N) = DF +D*(FPP(I) +B +B)/6. 

IND = 1 

IF <M0DE) 81,81,61 
FPPP(J) = FQM 

V = FPP(J) 

I = J 


J = J +K 

DS = S(J) -S(I) 

U = FPP(a) 

FPPP(J) = FPPP(I) +.5ADS*(FU) +F(a> 

V = U 

IF (J -N) 71,81,71 

81 RETURN 
END 


+U)/6. 


-DSADSA(U 


+V)/12.) 


SUBROUTINE INTPL(MI,NI,SI,FI,M,N,S,F,FP,FPP,FPPP,HODE) 

C IMPLICIT REALA8 (A-H,0-Z) 

INTERPOLATION USING TAYLOR SERIES 

ADDS CORRECTION FOR PIECEWISE CONSTANT FOURTH DERIVIATIVE 
IF MODE GREATER THAN 0 

DIMENSION SI(1>,FI(1),S(1),F(1),FP.(1),FPP(1),FPPP(1) 

K = IABS(N -M) 

K = (N -M)/K 

I = M 

MIN = MI 

NIN = NI 

D = S(N) -S(M) 

IF (DA(SKNI) -SKMI))) 11,13,13 
11 MIN = NI 

NIN - HI 

13 KI = IABS(NIN -MIN) 

IF (KI) 21,21,15 
15 KI = (NIN -MIN)/KI 

21 II = MIN -KI 

C =0. 

IF (MODE) 31,31,23 
23 C =1. 

31 II = II +KI 

SS = SKID 

33 I = I +K 
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O O o 



IF (I 

-N) 35,37,35 

IF (DA(S(I) -SS)) 33,33,37 

J 

= I 

I 

= I -K 

SS 

= SS -S(I) 

FPPPP 

= C*(FPPP(J) -FPPP(I))/( 

FF 

= FPPP(I) +.25ASSAFPPPP 

FF 

= FPP(I) +SSAFF/3. 

EF 

= FP(I) +.5ASSAFE 

FKII) 

= F(I) +SSAFF 

IF (II 

-NIN) 31,41,31 


41 RETURN 
END 



TURBINE GRID (PROFILE) 
FSYM ENU FNL 


0.0 

19. 

19. 

TOTO 

COVO 

PIO 

32. 

0. 

0.000 

FHITl 

FMIT2 

HUAL 


1.0 

1.0 

ALl 

AL2 

0.000 

0. 

TRAIL 

SLOPT 

5.0 

0.0 

XP(I) 

YP( I) 

0.000000 

0.000000 

0.015700 

0.027706 

0.108800 

0.088597 

0.212100 

0.137853 

0.315400 

0.180847 

0.399900 

0.212482 

0.496500 

0.242702 

0.568900 

0.259087 

0.629300 

0.266360 

0.689700 

0.264002 

0.738000 

0.250662 

0.786200 

0.222683 

0.822100 

0.192298 

0.849100 

0.167692 

0.885000 

0.133698 

0.911800 

0.108015 

0.944000 

0.079050 

0.976200 

0.053698 

0.999000 

-0.002591 

XP(I) 

YP(I) 

0.000000 

0.000000 

0.009900 

-0.044558 

0.037800 

-0.059870 

0.123200 

-0.072826 

0.208700 

-0.082956 

0.294100 

-0.089020 

0.379500 

-0.090128 

0.447900 

-0.086748 

0.533300 

-0.077099 

0.601600 

-0.064913 

0.670000 

-0.049004 

0.738300 

-0.030830 

0.789500 

-0.016857 

0.840800 

-0.004127 

0.875000 

0.001918 

0.926200 

0.003198 

0.960300 

-0.006010 

0.974700 

-0.012163 

0.999000 

-0.002591 

PRl 

PR2 

1. 

1.0 


0.220 

STEP 

0 . 

XSING 

.001050 


BUNCH 

.0 


FNX 

52. 

P20 

0.00 


FMY 

IG. 

P30 

0.0 


FMESH 

1 . 

P40 

0.00 


VIS 

1.0 

P50 

0.00 


FMl FM2 DHACH 

.500 0. 0. 

YSING 

0.000 


IRD lURI IGSR (3F10.6,3I&) 
-2 -1 +0 


QC 

0.90 

GHESH 

1.0 


22 . 


o o o o o 


Appendix II 


PARAMETER(KX =130,KY=67) 

C PROGRAM GRID2 (FOR CASCADE GEOMETRIES-USES OUTPUT FROM GRIDl 

DIMENSION XVM(KX,KY),YVM(KX,KY),DXP(KX),DYP(KX),C(KY) ,FP(KY),FM(KY) 
DIMENSION XVP(KX,KY),YVP(KX,KY),XN(KX,KY),YN(KX,KY) 

DIMENSION XF(KX,KY),YF(KX,KY),FI(KX),FJC(KX),FIC(KX) 

DIMENSION XU(KX),YCU(KX),YCL(KX) 

DIMENSION XI(10>,YI(10),XII(5>,YII(5),El(10),E2a0),E3(10) 

DIMENSION YVPN(KX,KY) ,YNN(KX,KY),F1(KX) 

DIMENSION XFF(KX,KY),YFF(KX,KY) 

DIMENSION Y01(KX),Y02(KX),YU(KX),YL(KX) 

PI=4.*ATAN(1,0) 
lURIT = 1 
DELTG = 0.0461 
8SCAL = 2.53 
DIR = 180. /PI 
THETA = 32.25/DTR 
ALPHA=-2.5/DTR 
A1 = 1.0113 

A2 = 0.0000 

A3 = -0.4240 

A4 = -0.4484 

C THE NEXT FOUR LINES DEFINE LOCATION OF THE CASCADE 

JBODY = 9 
JBP = JBODY +1 
ILE = 15 
ITE = 39 

IR = 1 
NPLT= 1 
STR = 0.5 

INPUT PROGRAM INITIALISATION PARAMETERS 
IREV =1 STATOR GRID, 0 -ROTOR GRID 

NAV = NO. OF TIMES SMOOTHING APPLIED (ABOUT 25 TO 50) 

ICR = 0 COARSE GRID 1 FINE GRID (UPSTREAM) 

INPUT FROM UNIT 8 THE TWO DIMENSION GRID GENERATED BY GRIDl. FOR 
READ(5,J0 IREU,NAV,ICR 
READ(8,A) NX, NY 
DO 20 J=1,NY 
DO 20 1=1, NX 

READ(8,A)XVP(I,3),YVP(I,J) 

XVP(I,J)=XVP(I,J)/GSCAL 
YVP(I,J)=YVP(I,J)/GSCAL 
20 CONTINUE 

C GENERATE THE CARTESIAN GRID WITH CORRECT OUTER PERIODIC SHAPE 

C THIS IS USED TO BLEND WITH THE GRID INPUT EARLIER 

XB6 =XVP(1,1) 

YBG =YVP(1,1) 

DELTX =ABS(XVP(1,1)-XVP(2,1))A.93 

DO 25 1=1, NX 

YN(I,1) = YBG 

YNN(I,1) = YBG 

XN(I,1) = XBG +DELTXA(I-1) 

25 CONTINUE 

DO 30 J=2,NY 

DELTY =-ABS(YVP(l,l)-YVP(l,2)) 
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JJ =J 
Cl = 1.38 
C1R= 1.1 

IE (J.GE.JBP) THEN 
Cl = .82 
CIR =1.1 
END IE 

IE (J.EQ.JBP) DELTY = 0.00000000 
IE (J.GE.JBP) JJ=J-1 
DO 30 1=1, NX 
XN(I,J)=XBG+DELTXA(I-1) 

YN< I,J)=YN( I,J-1) +DELTYA0.773AC1 
YNN( I, J)=YNN( I, J-1) +DELTYA0.773AC1R 
30 CONTINUE 

C DEEINE BLENDING EUNCTIONS IN I AND J DIRECTIONS 

ITE1C=ILE 
ITE1S= ILE/2 +1 
ITE2C=ITE 
DO 35 I=1,ITE1C 

EIC(I)=0.5A(1.-C0S(EL0AT(I-1)API/EL0AT( ITEIC-D) ) 

El< I) =0.5A(1 .-COS( FLOAT ( I-ITE1S)API/EL0AT( ITEIC-ITEIS) ) ) 
IE (I.LE.ITEIS) El(I) =0.0 
35 CONTINUE 

DO 95 J=1,NYC2 

EJC(J>=0.5A(1.-C0S(FL0AT(J-1)API/EL0AT(NYC2-1) )) 

95 CONTINUE 

DO 98 J=NYC2+1,NYC 
J1=NYC+1-J 

EJC(J)=0.5A(1.-C0S(EL0AT(J1-1)API/EL0AT(NYC-NYC2-1))) 

C FJC<J)=1. 

98 CONTINUE 

ITE2M = ITE2C 
DO 45 I=ITE1C,ITE2M 
Fl(I) = 1. 

EIC(I)=1. 

45 CONTINUE 

ITP = ITE+I 

DO 80 I=ITE2M+1,NX 

I1=NX-I+1 

EIC(I)=0.5A(1.-C0S(EL0AT(II-1)API/EL0AT(NX-ITE2M))) 

FKI) = EIC(I) 

80 CONTINUE 

C THE SPLITTING OE GRID LINE AHEAD OE LEADING EDGE AND 

C BEYOND TRAILING EDGE 

FACT = 0.5 

IE (ICR.EQ.l) FACT =.25 
DELTYG = -DELTGAEACT 
IIE2M = ITE2C-1 
IFF = ITP 
XILE = XVP(ITEIC,JBODY) 

YILE = YUPaTElC,JBODY) 

DO 85 I=ITE2M,IEE 
I1=ITP-I+1 

EI(I)=0.5A(1.-C0S(FL0AT(I1-1)API/EL0AT(ITP-ITE2M))) 
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85 CONTINUE 

DO 88 J=1,NY 
lEACT = -1 

IF (J.GE.JBP) IFACT = 1 

XVPdIElCjJ) = 0.5*(XVP(IIE1C,J)+XVP<ITE1C+1,J)) 

DO 11 I=1,ITE1C 

FACT = FIC(I)A0.75 

IF (IREV.EQ.O) FACT =F1(I)A0.75 

YVP(I,J) = YVP(I,J)+DELTYGAIFACTA(1.-FACT) 

YNN(I,J) = YNN(I,J)+DELTYGAIFACTA(1.-FACT) 

YN(I,J) = YN(I,a)+DELTYGAIFACTA(l.-EACT) 

11 CONTINUE 
lEND = ITP 

IF (IREV.EQ.O) lEND =NX 
DO 12 I=ITE2C-1,IEND 
FACT = FKI) 

YVP(I,J) = YVP(I,J)+DELTYGAIFACTA(1.-FACTA0.75) 
YNN(I,J) = YNN(I,J)+DELTYGAIFACTA(1.-FI( DA0.75) 

12 YN(I,J) = YN(I,J)+DELIYGAIFACIA(1.-FI(I)A0.75) 

88 CONTINUE 

C INTERPOLATE FOR BODY POINTS ON EITHER SIDE OF L.E AFTER 

C SPLITTING 

XKl) = XVP(ITE1C+2,JB0DY) 

YKl) = YVP(ITE1C+2,JB0DY) 

XK2) = XVP(ITE1C+1,JB0DY) 

YK2) = YVP(ITE1C+1,JB0DY) 

XK3) = XILE 
YK3) = YILE 
XK4) = XVP(ITE1C+1,JBP) 

YK4) = YVP(IIE1C+1,JBP) 

XK5) = XVP(ITE1C+2,JBP) 

YI(5> = YVP(ITE1C+2,JBP) 

YII(1)= YVP(ITE1C,JB0DY) 

YII(2)= YVP(ITE1C,JBP) 

MI =1 

CALL SPLIF(1,5,YI,XI,E1,E2,E3,2,0.,2,0.,0,0.,IND) 

CALL INTPL(MI,2,YII,XII,1,5,YI,XI,E1,E2,E3,0) 
XVPdTElC.JBODY) =XII(1) 

XVPdTElC,JBP)=XII(2) 

143 CONTINUE 

DO 90 J=1,NY 
DO 90 1=1, NX 

YNMd,J)=YNNd,J)+DELTYA1.98 
90 CONTINUE 
NYC2=NY/2 
NYC=NY 

C GRID BLENDING 
DO 100 3 =1,NY 
DO 100 I =1,NX 

FU = FJC(J)AFJC(J)AFICd)AFICd) 

FU = FUAFU 

XN(I,J) = (l.-FU)AXNd,J)+FUAXVP(I,J) 

YNd,J) = (l.-FU)AYNd,J) + FUAYVP(I,J) 

YNN<I,3) = (l.-FU)AYNNd,3)+FUAYVP(I,J) 
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100 CONTINUE 

DO 105 J = 1,NY 

IF (J.GT.JBODY) THEN 

JM = NY 

JB =JBP 

ELSE 

JM =1 

JB =JBODY 

END IF 

VAL = 0.4 

DO 105 1=1, NX 

CON = (YN(I,JM)AVAL-YN(I,JB))/(YN(I,JM)-YN(I,JB)) 
YVP(I,J)=YN(I,JB)+CON*(YN( I,J)-YN(I,JB)) 

CON = (YNN(I,JM)*VAL-YNN(I,JB))/<YNN(I,JM)-YNN(I,JB)) 
YVPN(I,J)=YNN(I,JB)+CONA(YNN(I,J)-YNN(I,JB)) 

105 CONTINUE 

DO 110 J=1,NY 
DO 110 1=1, NX 
YN(I,J)=YVP(I,J) 

YNN(I,J)=YVPN(I,J) 
no CONTINUE 

IE (IREV.GE.5) GO TO 2223 

C ADD CAMBER LINE DESCRIBED BY Y = A1+A2AX+A3AXAX+A4AXAXAX 
DO 115 J=1,NY 
DO 115 1=1, NX 
XU(I)=XN(I,J) 

YCU( I)=A1-A3AXU( I)AXU( I)-A4AXU( I)AXU( I)AXU( I) 
YVP(I,J)=YNN(I,J)+YCU(I) 

YN(I,J)=YN(I,J)+YCU(I) 

YVP(I,J)=-YVP(I,J) 

XVP<I,J)=XN(I,J) 

XVM(I,J)=XVP(I,J) 

YVM(I,J)=YVP(I,J) 

115 CONTINUE 

DO 120 J=1,NY 
JJ = NY+l-J 
DO 120 1=1, NX 

XVP(I,JJ)=XVM(I,J) 

YVPCI,JJ)=YUM(I,J) 

120 CONTINUE 

C ROTATE THE ROTOR GRID BY THETA ABOUT L.E 

XU = (XVP(ITE1S,JB0DY>+X0P( IIE1S,JBP))A0.5 

Yll = <YVP(ITE1S,JBODY)+YVP(ITE1S,JBP))A0.5 

DO 125 J=1,NY 

JB = JBODY 

IE (J.GE.JBP)JB =JBP 

DO 125 1=1, NX 

XI =XVP(ITE1S,J) 

Y1 =YVP(ITE1S,J) 

IF <J.NE.JB> GO TO 222 
IF <I.GE.ITE1C.AND.I.LE.ITE2C) THEN 
XI = Xll 
Y1 = Yll 
END IF 
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222 


CONTINUE 
XV =XVP(I,a) 

XVP(I,J)=X1+(XVP(I,J)-X1)*C0S(THETA)+(YVP(I,J)-Y1)ASIN(THETA) 
YVP(I,J)=Y1+(YVP(I,J)-Y1)*C0S(THEIA)-(XV-XI)ASIN(THETA) 

125 CONTINUE 

C GRID SMOOTHING 

DO 128 lAV =1,NAV 
DO 121 3=1, NY 
DO 121 1=1, NX 
XVM(I,J)=XN(I,J) 

YVM<I,J)=YN<I,J) 

121 CONTINUE 

3BMH = JBODY -1 

JBPP = JBP+1 

JNORM= NY -JBODY 

DO 122 3=2, NY-1 

IF (3. LE. JBODY) THEN 

ALPHA =1. -FLOAT (JBMM-J)/ELOAI(JNORM) 

ELSE 

ALPHA =1.- ELOAT(J-JBPP)/FLOAT(JNORM) 

END IF 

IF (J.EQ.JBODY.OR.J.EQ.JBP) GO TO 122 
DO 122 1=2, ITP 

XN(I,J)=(XVM(I+1,J)+XVM(I-1,J)+ 

1ALPHAAXV«(I,J+1)+XVM( 1, 3-1 )*ALPHA)/(2.+2.AALPHA) 
YN(I,3)=(YVM(I,J+1)+YVM(I,J-1))*.5 

122 CONTINUE 

126 CONTINUE 
ILEl =ILE -1 
3BO H JBODY 

XAV = 0.5A(XN(ILE1,JB0)+XN(ILE1,JB0+1)) 

YAV = 0.5A(YN(ILE1,JB0)+YN( ILEl,JBO+D) 

XN(ILEl+l,JBO+l) = 0.25*(XN<ILE1+1, JB0+1)A2.+XN( ILEl, JBO+1) 
1+XN<ILE1+2,3B0+1)) 

YN(ILEl+l,JBO+l) = 0.25A(YN(ILE1+1,JB0+1)A2.+YN(ILE1, JBO+1) 
1+YN(ILE1+2,JB0+D) 

XAV = 0.5*(XVP(ILE1,JBO)+XVP(ILE1,JBO+1)) 

YAV = 0.5A(YVP(ILE1,JB0)+YVP(ILE1,JB0+1)) 

XVP(ILE1+1,3B0) = 0.25A(XVP(ILE1+1,JB0)A2.+XVP(ILE1,JB0) 
1+XVP<ILE1+2,3B0)) 

YVP(ILEl+l,JBO) = 0.25A(YVP<ILE1+1,3B0)A2.+YVP(ILE1,JB0) 
1+YVP(ILE1+2,JB0)) 

DO 135 lAV =1,NAV 
DO HI 3=1, NY 
DO HI 1=1, NX 
XVMa,3)=XVP(I,J) 

YVH<I,3)=YVP(I,3) 

HI CONTINUE 

DO 142 3=2, NY-1 
IF (3. LE. JBODY) THEN 
JBOO = JBODY 

ALPHA =l.-(FLOAT(JBMM-J)/JNORM) 

ELSE 

JBOD = JBP 
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ALPHA =1.-(FL0AT(J-JBPP)/JN0RH) 

END IF 

IF (J.EQ.JBODY.OR.J.EQ.JBP) GO TO 142 
DO 142 1=2, IIP 

XAV = (XVM(I+l,J>+XVM(I-l,J)+XVH(I,J+l)+XV‘<a,J-l))*0.25 
YAV = (YVM(I+l,a)+YVH(I-l,J)+YVM(I,J+l)+YVM(I,J-l))A0.25 
IF (lAV.GT.lO) ALPHA =0.0 
XVP(I,a)=ALPHAAXAV+(l.-ALPHA)*XVM(I,J) 
YVP(I,J)=(YVH(I,J+1)+YVM(I,J-1))*.5 
142 CONTINUE 

135 CONTINUE 

C STRETCHING OF STATOR GRID 

XTE =(XN(ITE,JB0DY)+XN(ITE,JBP))A0.5 

YTE =(YN(ITE,JB0DY>+YN(ITE,JBP))A0.5 

XN(ITE,JBODY) =XTE 

YN(ITE,JBOIiY) =YTE 

XN(IIE,aBP)=XTE 

YN(ITE,JBP)=YTE 

XTE =(XVP(ITE,JB0DY)+X0P(ITE,JBP))A0.5 

YTE =(YVP(ITE,JB0DY)+YVP{ITE,JBP))A0.5 

XVP( ITE,JBODY) =XTE 

YVP(ITE,JBODY) =YTE 

XVP<ITE,JBP)=XTE 

YVP(ITE,JBP)=YTE 

ILM = ILE-7 

DO 300 3=1, NY 

DO 300 1=1, ILE 

DELTI = FLOATULE-D/FLOAKILE-U 
XN(I,J) = -STR*DELTIADELTI*DELTI+XN(I,J) 

300 CONTINUE 

DO 16 3=1, NY 

SL0P=(YVP(ITE2C-2,3)-YVP(ITE2C-l,3))/ 

1(XVP( ITE2C-2,3)-XVP( ITE2C-1,3) > 
YVP(ITE2C,3)=YVP(ITE2C-1,3)+SL0PA(XVP(ITE2C,3)-XVP(ITE2C-1,3)) 
SL0P=<YN(ITE2C-2,3)-YN( ITE2C-1 ,3) )/ ( XN ( ITE2C-2, 3)-XN( ITE2C-1 ,3) ) 
16 YN(IIE2C,3)=YN(ITE2C-1,3)+SL0PA(XN( ITE2C, 3)-XN(ITE2C-l ,3) ) 

C SHIFT THE ROTOR GRID TO GET COMBINED GRID 

IF (IREV.EQ.l) GO TO 2223 
XVC = XVP(ILM,1) 

XC = XN(ITP,1) 

XDIF = XC - XVC 
IF (IREO.EQ.O) GO TO 2221 
DO 170 3=1, NY 
YC = YN(ITP,3) 

YVC = YVP(ILM,3) 

DO 170 I = ILM, NX 
XVP(I,3) = XVP(I,3)+XC-XVC 
170 CONTINUE 

DO 175 3 =1,NY 
XN(ITP,3)=XN(ITP,1) 

XC = XN(ITP,3) 

YC = YN(ITP,3) 

XVC = XVP(ILM,3) 

YVC = YVP<ILM,3) 


DO 175 I=ILM,NX 
YVP(I,J)=YVP(I,J)+(YC-YVC)+DY 
175 CONTINUE 

NXN = ITP+l+NX-ILM 
DO 333 J=1,NY 
11=1 

DO 333 I=ITP+1,NXN 
IN =ILH+II 
XN(I,3)=XVP(IN,J) 

YN(I,3)=YVP(IN,J)+DY2 
II =11+1 
333 CONTINUE 

NX =NXN-1 
GO TO 2223 
2221 CONTINUE 

DO 2222 J=1,NY 
DO 2222 1=1, NX 
II = 1+7 

XN(I,J)=XVP(I1,J)+XDIE 

YN(I,J)=YVP(I1,J) 

2222 CONTINUE 

2223 CONTINUE 

IF (IREV.LE.l) THEN 
NX =ITP 
END IF 

C INTERPOLATE (LINEAR) TO FINE GRID FOR ICR 

IF (ICR.NE.l) GO TO 610 
KYU = JBODY 
KYB = JBP 
JJ = 0 

DO 500 3 =1,KYU 
33 = 2*3-1 

3P = 33+1 
DO 500 1=1, NX 
XVP(I,J3) = XN(I,J) 

YVP(I,33) = YN(I,3) 

IF (J.EQ.KYU) GO TO 500 
XVP(I,3P) = (XN(I,3)+XN(I,3+1))*0.5 

YVP(I,JP) = 0.5*(YN(I,3)+YN(I,J+D) 

500 CONTINUE 

DO 600 3 =KYB,NY 
33 = 2*3-2 

3P = 33+1 
DO 600 1=1, NX 
XVP(I,33) = XN(I,3> 

YVP(I,33) = YN(I,3) 

IF (3.EQ.NY) GO TO 600 

XVP(I,3P) = (XN(I,3)+XN(I,3+1))*0.5 

YVP(I,3P) = 0.5*(YN(I,3)+YN(I,3+D) 

600 CONTINUE 

NY = 2*NY-2 
DO 610 3=1, NY 
DO 610 1=1, NX 
XN(I,3) = XVP(I,3> 
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620 

C 


1200 


1250 

1275 


980 


YN(I,J) = YVP(I,J) 

CONTINUE 

IE (IREV.EQ.l) 60 TO 620 
DO 620 J=1,NY 
XN(NX,J) = XN<NX,1) 

XN(1,J) = XN(1,1) 

CONTINUE 

ADDITOINAL SMOOTHING FOR ROTOR GRID 
ISKIP = 1 

IF (IREV.LT.l) ISKIP = 0 
IF (ISKIP. NE.l) GO TO 1275 
JUP = JBODY 

IF (ICR.GE.l) JUP =2*JB0DY-1 
NXF = NX 

IF(IREV.GT.l) NXF = ITP 
DO 1200 1=1, NXF 
YT = YN(I,1)-YN(I,NY) 

YBl = YN(I,JUP) 

YB2 = YN(I,JUP+1) 

YOKI) = YN(I,1) 

Y02(I) = YN(I,NY) 

YU(I) = (YT+YBl+YB2>/2. 

YL(I) = (YBl+YB2-YT)/2. 

DO 1250 J=1,NY 

JBOD = JBODY 

IF (J.GE.JBP) JBOD = JBP 

IF (ICR.EQ.l) THEN 

JBOD = 2AJB00Y-1 

IF (J.GE.JBOD+1) JBOD = JBOD+1 

END IF 

DO 1250 1=1, NXF 
YBODY = YN(I,JBOD) 

YMUL = (YU<I)-Y01(I))/(Y01(I)-YB0DY) 

IF (J.GT.JBODY+1) YMUL =(YL( I)-Y02( I) >/(Y02( I)-YBODY) 

YN(I,J)= YN(I,J)+(YN(I,J)-YBODY)AYMUL 

CONTINUE 

CONTINUE 

JBO = JBODY 

IF (ICR.EQ.l) JBO = 2AJB0DY -1 
KY2 = NY/2 
IF (IREV.NE.O) GO TO 881 
ILR =8 

DO 881 N =l,NAV-40 
DO 980 J=1,NY 
DO 980 1=1, NX 
XVM(I,J) =XN(I,J) 

YMM(I,J) =YN(I,J) 

CONTINUE 
DO 880 J=2,NY-1 
JP = 1 
JBOD = KY2 

IF (J.GT.KY2) JBOD = KY2+1 
IF (J.GT.KY2) JP =NY 
DO 880 I=2,ILR-1 
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OMEGA = ABS(I-2)AAEiS(J-JP)/(ABS(2.-ILR)AABS(JP-JB0D)) 
XAV = 0.25A(XVH(I+l,a-l)+XVM(I+l,J+l)+XVM(I-l,J-l)+ 
1XVM(I-1,J+D) 

YAV = 0.25A(YVM(I,J-1)+YVM(I,J+1)+YVM(I,J+1) 

XN(I,J) = OMEGAAXAV+(l.-OMEGA)AXN(I,a) 

YN(I,J) =OMEQAAYAV+(l.-OMEGA)*YN(I,J) 

880 CONTINUE 

881 CONTINUE 

C INTERPOLATE FOR 3 DIRECTION AND STORE GRID COORDINATES 

IE (IREV.GE.2.0R.IURIT.EQ.0) GO TO 750 
IF <IREV.EQ.l) THEN 
JREV = 15 

URITE (12,1004) NX, NY 
ELSE 

URITE (12,1003) NX, NY 
JREV = 8 
END IF 

YDIF = -YN(JREV,JBO) 

DO 2224 J=1,NY 
DO 2224 1=1, NX 
YN(I,J) = YN(I,J)+YDIF 
2224 CONTINUE 

RAD = 4.777 
RADI = 4.503 
DO 790 J =1,NY 
DO 790 I =1,NX 
790 YNN(I,J) = YN(I,J)/RAD 
DO 810 K = 1,15 
RAD = RADI 4(K-1)A0.0G8G 
URITE(12,1005) RAD,K 
DO 700 J=1,NY 
DO 700 1=1, NX 
X = XN(I,J) 

Y = RADASIN(YNN(I,J)) 

Z = RADACOS(YNN( I,J)) 

WRITE(12,1001) I,J,X,Y,Z 
700 CONTINUE 

810 CONTINUE 

750 CONTINUE 

1001 FQRMAT(2I8,8F12.6) 

1002 F0RMAT(8F12.6) 

1003 F0RMAT(40X,10HR0T0R GRID,5H NX =,I5,5H NY =,I5) 

1004 FORMAT(40X,11HSTATOR GRID,5H NX =,I5,5H NY =,I5) 

1005 F0RMAT(//40X,8HRADIUS =,F12.6,4H K =,I5) 

10 F0RKAT(2I8) 

40 F0RMAT(2F12.6,2I8) 

STOP 
END 

SUBROUTINE INTPL(MI,NI,SI,FI,M,N,S,F,FP,FPP,FPPP,MODE) 
INTERPOLATION USING TAYLOR SERIES 

ADDS CORRECTION FOR PIECEUISE CONSTANT FOURTH DERIUIATIOE 
IF MODE GREATER THAN 0 

DIMENSION SI(1),FI(1),S(1),E(1),FP(1),FPP(1),FPPP(1) 
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14 

U 

= 

-1. 


V 

= 

-DS*VM 


GO TO 25 



21 

I 

= 

J 


J 

= 

J +K 


DS 

= 

S(J) -S(I) 


IF (D*DS)81 

,81,23 

23 

DF 


(F(J) -F(I))/DS 


B 

= 

l./(DS +DS +U) 


U 


BADS 


V 

= 

BA(6.ADF -V) 

25 

FP(I) 

= 

U 


F?P(I) 

r 

V 


U 

= 

(2. -U)ADS 


V 

= 

6.ADF +DSAV 


IF (J -N) 21,31,21 

31 

IF (KN -2) 

32,33,34 

32 

V 

s 

(6.AVN -V)/U 


GO TO 35 



33 

V 

= 

VN 


GO TO 35 



34 

V 

=: 

(DSAVN +FPP(I))/(1 

35 

B 

= 

V 


D 

= 

DS 

41 

DS 

s 

S(J) -S(I) 


U 

s 

FPP(I) -FP(I)AV 


FPPP(I) 

s 

(V -U)/DS 


FPP(I) 

= 

U 


FPU) 

= 

<F<J) -F(I))/DS - 


V 

= 

U 


J 

= 

I 


I 

s 

I -K 


IF (J -M) 

41,51,41 

51 

I 

s: 

N -K 


FPPP(N) 

= 

FPPP(I) 


FPP(N) 

= 

B 


FP(N) 

= 

DF +DA(FPP(I) +B 


IND 

= 

1 


IF (H0DE)81, 

,81,61 

&1 

FPPP(J) 

= 

EQM 


V 

= 

FPP(J) 

71 

I 

= 

J 


J 

= 

3 +K 


DS 

s 

S(J) -S(I) 


U 

= 

FPP(J) 


FPPP(J) 

= 

FPPP(I) +.5ADSA(FC 


V 

= 

U 


IF (J -N) 


71,81,71 

81 

RETURN 




END 




+EP(D) 


-DS*(V +U +U)/&. 


+B)/6. 


+F(J) -DS*DS*(U 


+V)/12.) 


